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We present a collision model for particle-particle and particle-wall interactions in interface- 
resolved simulations of particle-laden flows. Three types of inter-particle interactions are 
taken into account: (1) long- and (2) short-range hydrodynamic interactions, and (3) solid- 
solid contact. Long-range interactions are incorporated through an efficient and second- 
order accurate immersed boundary method (IBM). Short-range interactions are also partly 
reproduced by the IBM. However, since the IBM uses a fixed-grid, a lubrication model is 
needed for an inter-particle gap width smaller than the grid spacing. The lubrication model is 
based on asymptotic expansions of analytical solutions for canonical lubrication interactions 
between spheres in the Stokes regime. Roughness effects are incorporated by making the 
lubrication correction independent of the gap width for gap widths smaller than ~ 1% of 
the particle radius. This correction is applied until the particles reach solid-solid contact. 

To model solid-solid contact we use a variant of a linear soft-sphere collision model capable 
of stretching the collision time. This choice is computationally attractive because it allows 
to reduce the number of time steps required for integrating the collision force accurately 
and is physically realistic, provided that the prescribed collision time is much smaller than 
the characteristic timescale of particle motion. We verified the numerical implementation of 
our collision model and validated it against several benchmark cases for immersed head-on 
particle-wall and particle-particle collisions, and oblique particle-wall collisions. The results 
show good agreement with experimental data. 
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I. INTRODUCTION 


Flows laden with solid particles appear widely in both nature and industry. Examples are the 
transport of sediments in a river, the enhanced mixing due to the presence of particles in a fluidized 
bed reactor, and the flocculation/sedimentation processes in the treatment of drinking water. In 
many cases the flow is turbulent, the size of the particles is comparable to or larger than the 
Kolmogorov length-scale (i.e., the particles have a finite-size), and the volume fraction of particles 
may be considerably high such that inter-particle interactions are dynamically important [l|. 

Studying flows laden with finite-size particles using interface-resolved direct numerical simula¬ 
tions (DNS) has recently become possible with the development of efficient numerical methods, 
such as immersed boundary methods (IBM) [^, together with the continuous increase in comput¬ 
ing power. Such simulations provide detailed insight in the flow dynamics at the particle scale 
and beyond. The governing equations for the fluid phase and the particles are directly coupled 
with each other through the no-slip/no-penetration condition at the particles’ surfaces (i.e., 2-way 
coupling) without the need of parameterizing the drag force between the phases. Also, long-range 
hydrodynamic interactions between particles (i.e., 4-way coupling) are naturally reproduced by 
these methods. However, when the particle volume fraction is high, additional models are re¬ 
quired to account for short-range hydrodynamic solid-solid interactions (lubrication forces) and 
solid-solid contacts. Otherwise, the realism of the simulation may be compromised by a poor de¬ 
scription of these interactions. For instance, by under-predicting lubrication-enhanced clustering 
of inertial particles, as observed for homogeneous isotropic turbulent flows Q]. The challenge is 
to find a model able to reproduce short-range particle-particle and particle-wall interactions with 
the required realism and with little effect on the computational efficiency of the overall numerical 
algorithm. 

We consider non-Brownian spherical particles, which are sufficiently large such that inter-surface 
forces as the Van der Waals force and the electrostatic double-layer force can be neglected [ 4 ]. Also, 
cohesive forces, which are relevant for wet granular media a. »e We ..e 

applicability of the model to cases where 4-way coupling is required, but where the solid volume 
fraction is not extremely high such that good description of the macroscopic outcome of the collision 
(i.e., relative velocity prior to and after contact) is sufficient to model the suspension dynamics. 

Much work has been done in modeling of inter-particle (or particle-wall) collisions. Discrete 
element methods (DEM) have been successfully used to account for inter-particle collisions in 
simulations of gas-solid flows where hydrodynamic interactions between particles are negligible 
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(e.g., B, fl and B)- These collisions are often referred to as dry collisions. More recently, 
some studies used these same collision models for reproducing particle-particle and particle-wall 
interactions in viscous liquids, commonly referred to as wet collisions. In this case, fluid effects 


B- 


such as added mass, viscous dissipation and history forces become important 

The lubrication effects cannot be resolved by the overall numerical method (not without resort¬ 
ing to excessive grid refinement). This lack of spatial resolution can be circumvented by a closure 
model for lubrication interactions based on analytical solutions of these interactions in the Stokes 
regime (e.g., Q, Q, Q and Q). 

Many studies used variants of the soft-sphere collision model of Cundall and Strack B to com¬ 
pute the contact forces, because of its computational advantages for simulating dense suspensions 


when compared to hard-sphere models 


14l |. In the soft-sphere model, the normal force acting 


on the particle during a collision is computed from an equivalent linear spring-dashpot system in 
which the spring stiffness and dashpot coefficients are parameterized as function of the particle’s 
elastic properties. A limitation of this approach when applied to particle-laden flows is that the 
collision must be resolved with a time step that can be several orders of magnitude smaller than 
the time step of the Navier-Stokes solver for the fluid flow. This happens because the characteristic 
time scale of solid-solid contact is in general orders of magnitude smaller than the smallest time 
scale present in the flow [iB ]- [iB- However, it is possible to artificially stretch the collision time 
to a multiple of the time step with which the particle motion is integrated. In some studies this 
was done by decreasing the value of the spring stiffness and checking resulting the collision time in 
a trial and error procedure [l7]. This approach was avoided by others, who prescribed the desired 
collision time and computed the cor resp onding collision parameters by solving the equations of the 


n 


n. 


harmonic oscillator (e.g., [l]|, [1^, |l3l|b 

Experimental studies have shown that the fluid effects in the normal collision of a sphere onto 
a plane wall can be quantified by an effective normal coefficient of restitution, e^, defined as 
the ratio of the magnitudes of rebound and impact velocities. In particular, when experimental 
data of en/cn^d (where is defined in an analogous way as Cn but for a collision in a dry 
system) are plotted against the particle impact Stokes number, St = {l/9)ppUpDp/yi (where pp, 
Up, Dp and p are respectively the particle mass density, impact velocity, diameter and the fluid 
dynamic viscosity), the datasets for different fluids and particle types collapse in the same curve 


1 ^. This suggests that Cn, en,d and St are key parameters to describe a head-on wet collision. 


Hence, reproducing this scaling is an important test for any numerical method for resolving the 
flow conforming a particle combined with a collision model should pass. Several authors have 
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been able to reproduce it with different methodologies for resolving the particle-fluid interface, 


such as tensorial penalty methods 
and 


n 


,n 


Lagrange multiplier-based methods |l8l | or IBM 


, Q, 


2ll|). However, this benchmark experiment relies in a definition of impact and rebound 
velocities, which vary significantly in these references [^. Hence, if one solely resorts to this simple 
benchmark for validating the head-on collision model without careful comparison with experimental 
data, it can happen that the definitions of the impact and rebound velocities determined from the 
numerical simulation are not consistent with the measured quantities. 

The complexity of the problem increases when the collision is oblique. In this case, the relative 
motion between the contact surfaces has a tangential component. Two different kinds of motion 
can occur between the surfaces in contact; rolling and sliding. Rolling occurs when a point of 
contact has zero relative velocity with respect to the contact surface, otherwise sliding occurs. 
Moreover, when a particle flowing through a viscous liquid approaches a planar surface obliquely, 
it experiences not only lubrication forces due to the squeezing motion of the fluid through the 
gap, but also forces and torques due to relative translational and rotational shearing (see 2^ for a 
review). Finally, the frictional resistance of the contact surface in the presence of a viscous liquid 


can change abruptly due to piezoviscous effects when smooth particles collide obliquely 


M- 


To the best of the our knowledge, Kempe and Frohlich [1^ report the only collision model 
validated against experimental data of oblique particle-wall collisions in viscous liquids and against 
bouncing trajectories of particles colliding onto a planar surface in a viscous liquid. The latter 
benchmark validation is particularly interesting to reproduce because it does not rely on definitions 
of impact and rebound velocity. It therefore gives a finer indication of the success of the model 
to reproduce the canonical case of a particle-wall collision than reproducing data of enl&n,d vs St. 


Kempe and Frohlich 


121 ] computed the normal collision force from a non-linear spring-dashpot 


system. This was done so that the force-displacement relation agrees with Hertzian contact theory. 
The collision time was stretched by using a numerical procedure to solve the resulting equations 
of the non-linear spring-dashpot oscillator. For the tangential component, they developed a model 
based on the assumption that, throughout solid contact, a particle either rolls or slides, depending 
on the particles’ incidence angle. Although the approach of considering pure rolling for small 
incidence angles does not reproduce collisions with recoil of the contact point, their methodology 
can be easily adapted to account for it. Even though their model is able to reproduce normal 
and oblique collisions in viscous liquids with satisfactory realism, the fact that it needs an extra 
iterative procedure to deal with the non-linear spring when computing the normal force may 
deteriorate its computational performance for denser concentrations. Furthermore, the force law 
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for the tangential component of the collision force depends on the particles’ incidence angle, which 
is difficult to interpret, e.g., for cases in which geometrical constrains force sustained contact. 

In the present study we present a new model for wet particle-particle and particle-wall collisions 
in fully-resolved simulations of particle-laden flows. We show that a simple variant of a linear spring- 
dashpot model capable of stretching the collision time B, Q suffices for computing contact 
forces. The contact model can be seen as a linearized version of Hertz contact theory, and its 
choice is motivated by a separation between the time scales of solid-solid contact and particle 
motion. The advantage of using this model is that its parameters can be analytically determined 
from well-documented material properties and a desired collision time, which is computationally 
attractive. Moreover, it accounts for stick-slip effects at the contact point without requiring explicit 
definition of impact and rebound angles. Oblique collisions with recoil are explicitly accounted for 
by using a tangential coefficient of restitution et as input parameter. This contact force model is 
implemented in an efficient and second-order accurate IBM for particle laden flows developed in 


24l | and combined with a physical model for lubrication interactions and roughness effects. We 


found the experimental data used by Kempe and Prohlich 


121 ) to validate their approach to be a 


good set of canonical tests for which a physically realistic collision model should pass. We therefore 
validated our collision model against those distinct experimental cases which include the trajectory 


25l | and data oblique particle-wall collisions [23). 


of a sphere colliding onto a planar surface in a viscous liquid B , head-on particle-particle collisions 


This article is organized as follows. Section [TI] presents a brief overview of the physics of dry 
collisions of elastic spheres (III Ajl followed by the description of the methodology for computing 
contact force/torques ()IIB|1 . Section Hill addresses the effects of the interstitial fluid in a wet collision 
and our modeling strategy for lubrication interactions. The numerical implementation is addressed 
in Section HYl Section |V] explores the consequences of excessive and insufficient stretching of the 
collision time and presents the validation of the model against several benchmark experiments. 
Finally, in Section IVTI the conclusions and outlook are given. 


II. DRY COLLISIONS 

A. Physics 


When head-on inter-particle collisions take place in the absence of a viscous fluid, kinetic energy 
is dissipated exclusively due the contact mechanics. This energy loss can be described by a dry 
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coefficient of restitution, defined as the ratio of the relative rebound velocity to the relative 
impact velocity. The collision is referred to as oblique when the particles approach each other with 
an incidence angle just prior to contact 4>in and bounce with a rebound angle 4>out, as illustrated 
in Figure [TJ From these, it is convenient to define effective angles of incidence and rebound. 



FIG. 1. Schematic representation of an oblique inter-particle collision. For the sake of clarity we considered 
in this figure a reference frame moving with the light grey particle, which implies that the velocities sketched 
are relative velocities. Fn and Ft denote the normal and tangential component of the collision force. 


respectively as. 


^in = = tan{(j)in), and 


u. 


tn^n 


= e„,rftan(0o„t), 

with the normal dry coefficient of restitution Cn^d given by 

'^out.n 


^nA — 


Maw et al 


Uir. 


( 1 ) 

( 2 ) 

( 3 ) 


26( 1 explored this problem in detail. They used Hertzian contact theory to obtain the 


normal component of the collision force and velocity. Moreover, they assumed particles of the 
same material for which the contact area consists of stick and slip regions, and that slip could be 
modeled by a constant coefficient of sliding friction, fic- Their results show that three different 
types of impacts can occur, depending on the value of the following normalized incidence angle, 

( 4 ) 


Win — - 


2 - 

and the material- and geometry-dependent parameter, 

/ 1 \ 1-u 

X = 1 + 


2 - 


( 5 ) 
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where v is the Poisson’s ratio and K the normalized particle radius of gyration = 2/5 for a 
homogeneous solid sphere). Figure [2] shows ipout as a function of 'i/m as computed in their model, 
where ipout is the normalized rebound angle, defined analogously to V’m as. 


'4^ out 


2 1 - 1 / 


/Uc 2 - 1/ 


^out- 


( 6 ) 



FIG. 2. (Color online) Numerical solution of the model of Maw et al. 2^ for collisions between glass 


spheres, compared with experimental data of Foerster et al. 


271 and the model of Walton 281. The curve 
and rescaled to ipout vs tpin with 


and experimental data were extracted from a curve 4'o«t vs of 
the parameters of their homogeneous 3mm glass spheres, ly = 0.22 and = 0.092. The two vertical dotted 
lines delimit the three different types of impact and are given by tpin = 1 and tl’in = 4y — 1 Ri 4.2. 


The numerical solution of this model yields three distinct regions denoted in Figure[2]by I, II and 
III. First, for small incidence angles, < 1, the sphere sticks during contact because the normal 
component of the load is much larger than the tangential component. When the contact surface 
starts to shrink, small regions of micro-slip may occur due to tangential elastic recovery, which can 
spread throughout the entire contact area, leading to gross slip. Second, for an intermediate range 
of incidence angles, 1 < 'tpin < 4% — 1, the collision starts with gross slip, but the frictional stresses 
retard the tangential velocity, which rapidly drops to zero in the entire contact area (full stick). 
Finally, for higher incidence angles, 'ipin > 4% — 1, the tangential component of the load is even 
higher and gross slip occurs throughout the contact time. 


Walton 


28( 1 proposed a simplified hard-sphere model with three parameters: (1) a normal 


coefficient of restitution, (2) a tangential coefficient of restitution for non-sliding contact, 
defined as, 

_ '^out,t _ ^out . 




'^ind 




(7) 


and (3) a coefficient of sliding friction, ^c, to model the tangential component of the load. Ft, 
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acting on the particle when it is sliding, 

Ft = -fic\Fn\, ( 8 ) 

where Fn is the normal component of the contact force acting on the particle. This model assumes 
that the collision force acts at a single point and can be decomposed into a normal and tangential 
component. It further assumes that throughout the collision time the regime is either full stick or 
gross slip. From these three parameters, one can define the two lines which dictate the collision 
regime: 


= 


(stick), (9a) 

^m-Atc(l + l/-fii^)(l + e„,d),^m > (sHp), (9b) 

where Eq. (|9ap is obtained directly from the definition of and Eq. (I9bh by applying the 
definition of coefficient of sliding friction to relate the normal and tangential momentum impulses. 
Tt is the incidence angle above which the collision regime changes from full stick into gross slip: 


'ko.t('k*+) = ^ 1 + 


1 \ 1 “ 1 “ ^n,d 


The two models differ most significantly in the intermediate region of incidence angles, for which 
there may be periods of full stick and gross slip throughout the contact. Despite these differences, 
the s imp lified approach is able to reproduce experimental data reasonably well, as shown, e.g., in 
28l |. 27| (Eigure[2]) and 231]. The minimalistic nature Walton’s model makes it an attractive for 
problems where a detailed description of the contact mechanics is not required, which is in general 
the case for particle-laden flows. 


B. Modeling 

Legendre et al. [l^ demonstrated that collisions of spherical particles in a viscous liquid have 
a contact time larger but of the same order than the contact time in a dry system, predicted by 
Hertzian contact theory. They show that this contact time is four to five orders of magnitude 
smaller than the viscous relaxation time of the particle, depending on the impact Stokes number. 
This means that that the particle experiences a collision as a discontinuity in its motion. Even if 
the characteristic time scale of the particle motion is not dictated by the viscous relaxation time, 
(e.g., due to geometrical constrains in a flow with high volume fraction of particles) this clear 
separation of time scales typically remains valid. Hence, we require that the collision dynamics 
are realistically reproduced from a macroscale perspective, i.e., realistic approach and rebound 
velocities and timescale small enough such that this separation of time scales is satisfied. Hence, 
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it is convenient to use a model capable of stretching the collision time, so that that the overall 
numerical algorithm is not significantly penalized by the overhead introduced by the integration of 
the particles’ equations of motion. Furthermore, it is convenient to use a model with parameters 
that can be easily measured experimentally, such as the parameters of Walton’s model. Joseph and 
Hunt 2^ successfully used this model to describe experimental data from wet oblique collisions of 


spherical particles onto planar surfaces, which further supports its validity to describe collisions in 
a viscous liquid. 


I- 0 , 


81] to 


We found the variant of the soft-sphere contact model of Tsuji et al. 17|, described in 
be suitable for our problem. This approach has computational advantages for dense suspensions 
when compared to other alternatives such as hard-sphere models and allows the collision time to 
be stretched. The model consists on a linear spring-dashpot system in the normal and tangential 


directions, with a Coulomb friction slider in the latter, as sketched in Figure 3(a) In the following 
we describe the model, with differences in the definition of the tangential unit vector and in the 
value to which the tangential displacement is saturated. Figure [3(b) j illustrates the notation and 
reference frame adopted. 




FIG. 3. (a) Linear spring-dashpot model, (b) Notation and reference frame adopted for an inter-particle 
collision. 


The normal force acting on particle i due to a contact with particle j, with a relative velocity 
at the contact point given by 


Uij = (uj -b RiUi X Uij) - (uj -b RjiUj X Uji ), 


( 11 ) 


is the component of the collision force that acts along the direction of the line-of-centers (Figure 
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3(b)), 


riij = 


( 12 ) 


This collision force depends on the overlap distance between the two particles, 

^ij,n — {Ri T Rj ll^i ^*11) ) 

and the normal relative velocity of the contact point, 

^ij,n — (^jj ' 

and is obtained from the equivalent linear spring-dashpot system: 

^ij,n — kndij^ji 

where kn and rjn are the normal spring and dashpot coefficients, respectively. These are computed 
by solving for the motion of a linear harmonic oscillator Q] , and requiring that there is no overlap 
at the end of the collision, t = NAt, 


(13) 


(14) 


(15) 


■ ^ij') \t=NAt — 0, (f-b) 

and that the velocity at the end of the collision is given by the definition of en,d, 

{^ij,n ■ H-ij) \t=NAt — (^n,d {^ij,n ' |4=0' (1^) 

Note that we define the collision time, T„, as a multiple N of the time step of the overall numerical 
algorithm. At. This is convenient because - as our results will show - the outcome of a numerical 
simulation of a wet collision is more realistic if the fluid is allowed to adapt itself to the sudden 
changes in particle velocity. In practice, because r„ should be fixed during a collision, and At may 
vary in agreement with the stability criterion of the fluid solver, one should define the collision 
time as a multiple of the estimated time step of the numerical algorithm. 


The coefficients read. 


where 


me (vr^ + In^ en,d) 
(iVAt)2 


2me In en,d 
(NAt) 


me 

is the reduced mass of the particles. 



+ m. 


-1 


(18) 


(19) 


This approach can be seen as a linearized version of Hertzian contact theory. Since we model the 
collision as a discontinuity in the particle motion, it is sufficient to guarantee that the conditions 
specified in equations (I16h and (I17p are fulfilled and NAt is small enough so that the separation 
of time scales is satisfied in good approximation. One advantage of using a linear system is 
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that the spring and dash-pot constants can be determined analytically and a priori, which is 
computationally attractive. Notice that increasing the value of T„ reduces the spring stiffness, 
which makes the contact softer. This implies that excessive stretching of the collision time results 
in a large overlap between solid surfaces and consequently in a unrealistic delay of the particle 
rebound. On the other hand, the collision time should be sufficiently stretched so that the collision 
force is accurately resolved in time. We require that the maximum particle overlap, which is 
reached when the particles have zero relative velocity, = Sij^n\uij „=o^ is much smaller than 


Dr. 


T T* — o _ — _ 

\^ij * ^ij j |i=0 


(arcsin(7r/a)/7r) 


( 20 ) 


where a = + ln^(en,d) [81. Alternatively, if applicable, one can require that the maximum 

overlap due to the particle’s submerged weight = |1 — Pf/pp\g/kn) is much smaller than 


Dp: 


Tn « = 


I Dp a? 


9 \'^-pf/pp 


( 21 ) 


The tangential force is obtained analogously to but now with a Coulomb friction model 

to account for sliding motion: 


= min (II - hdij^t - II “ /^cF 


'‘V 


( 22 ) 


where 8ij t is the tangential displacement and tij the unit vector with the direction of the test 


force: 


tij — 




Hi. A m 11- (23) 

The coefficients kt and rjt are obtained in an analogous way by solving an harmonic oscillator for 
the tangential direction, and requiring that the definition of the tangential coefficient of restitution 
is fulfilled, 

(ujj^t ■ tjj) \t=NAt — (Wjj ■ tjj) |t=0) (2^) 

and that the collision times in the normal and tangential directions match (Tt = Tn). The values 
of the coefficients read, 

^ _ me,i (tt^- h In^ 2mejlnet,rf 

/ArA ,\0 1 Pt — t -\T \ J.\ ’ 


(NAt)^ 

where the reduced mass of the system is given by: 


(NAt) 


me,i = (l + l/AT^) ^ me. 


( 26 ) 


The tangential displacement of the contact point must be integrated in time from the imminence 
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of contact. From the integration of the relative tangential velocity at the point of contact we get 


S 


*nH-l 

ij,t 


fU+l 




u 


v,t 


dt, 


(27) 


where R is a rotation tensor which rotates ^ to the new local coordinate system at time level 
n + 1. 


The tangential force becomes independent of the tangential displacement of the spring when 
the particle starts sliding (Eq. (|22p l. If the tangential displacement is further incremented when 
the particle starts to s 
quently to sticking 


ide, unrealistic results can be obtained if the collision regime changes subse- 


29l |. Hence, the tangential displacement must be saturated in order to comply 


with Coulomb’s condition, whenever the collision is in the sliding regime 




xn+l _ 
^ij,t - 






(28a) 

(28b) 


After computing the contact forces acting at the point of contact, we determine the equivalent 
force and couple acting in the particle centroid: 

Fjj = + Fij^nj (29) 

= Rp (riy X Fij^t ). (30) 

The total collision force and torque are the sum of contributions of all the particles in direct contact 
with the particle i: 

= (31a) 

3 

Tf = ^Tf^.. (31b) 

3 

A wall is treated as a semi-inhnite spherical particle, which makes a particle-wall collision the 
limit case of a spherical particle with finite-radius, Ri, colliding onto a sphere with radius Rj ^ oo. 
Thus, the parameters for particle-wall collisions are computed in a similar way by taking this limit. 
The reduced mass is now given by me = and the normal overlap by Siw,n = (Ri — | |xi — 11) 

where rij^ is the unit-vector perpendicular to the wall and the coordinate of the point of contact 
on the planar surface. 



13 


III. EFFECTS OF THE INTERSTITIAL FLUID 

A. Lubrication effects 


A particle immersed in a viscous liquid experiences lubrication effects when moving close to and 
with finite relative velocity to another particle or wall. Assuming a drainage of the intervening 
liquid film in the Stokes regime, the force acting on the particle has an analytical solution that 


diverges when the non-dimensional gap-width, e = 5ij^n/Rp, tends to zero 3]|. Our IBM is able 


to reproduce this and other analytical solutions until a certain (small) value of e. For smaller 
gap-widths (< Ax) the IBM under-predicts this force due to a lack of spatial grid resolution. 
An approach that has been adopted for these cases is to keep the grid fixed and use lubrication 


models based on asymptotic expansions of analytical solutions for the 


ubrication force in the Stokes 


regime to compensate this lack of spatial grid resolution (e.g. [10|, [ll|], [1^ and [1^). Taking these 
effects into account has been proved to be important for computing realistic bouncing velocities in 
simulations of head-on particle-wall collisions in viscous fluids [l^ . 

Lubrication theory shows that ideally smooth particles would not reach actual solid-solid con¬ 
tact. Even if one accounts for the particles’ surface deformation due to the abrupt increase of the 



mental data for wet head-on collisions of a spherical particle onto a planar surface when a > h^- 
They argued that the contact occurs through the asperities, which are irregularly oriented, before 
elastohydrodynamic lubrication effects become important. This reasoning validates the approach 
used by several authors of setting the lubrication correction to zero for small gap-widths (e.g., [l^ . 
12 1 ) or making it independent of the gap-width {ll|. 


The most important component of the lubrication forces acting on the particle is the squeezing 
force acting along the line-of-centers, because its dominant term is oc 1/e in contrast to translational 
and rotational shearing, which diverge slower (oc Ine) and even for a value of e compliant with 
surface roughness have a negligible effect in the particle dynamics. Test simulations showed that 
the latter mentioned lubrication corrections had little effect on the results for immersed oblique, 
particle-wall collisions and therefore we decided to neglect them in the present study. 


We use a two-parameter model to account for normal lubrication interactions and roughness 
effects, as illustrated in Figure 01 When a spherical particle approaches a planar surface/another 
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particle, for a certain gap-width, the IBM cannot resolve the lubrication force acting on the 
particle. Hence, for gap-widths smaller than eax, we correct the lubrication force acting on the 
particle by adding to the Newton-Euler equations AFiub = —dn~ ^i^Ax)), where the 


Stokes amplification factor A is given by 


M- 



Lubrication fully resolved by IBM 


, IBM + 

/ Lubrication correction 


1BM + 

Lubrication correction 
with roughness effects 


Contact model takes over 


o 


FIG. 4. Schematic representation of the lubrication model. We illustrate the case of particle-wall interactions 
for the sake of simplicity. The model is analogous for particle-particle interactions. 


19 3 

Xvvi^) = -Ine-elne -b 0(1), 

^ 2e 20 56 ^ ^ 

Apu;(e) = - --Ine -—elne-b 0(1), 

£ t) Z L 


(32) 

(33) 


for lubrication interactions between two equal spheres, and between a sphere and a planar surface, 
respectively. 

The value o f ea x can be determined by simulating the slow approach of a sphere towards a 


planar surface 


3ll | or between two spheres 35|and determining up to which point the IBM is 


able to reproduce the lubrication interaction 


lll |. We illustrate this in Figure [5] by comparing 


the analytical solution with the simulations without lubrication correction and with lubrication 
correction. The corresponding values of eax for two different spatial resolutions are given in Table 
m To account for surface roughness, we saturate the Stokes amplification factor for gap-widths 


TABLE I. Parameters for the lubrication model 


Dp/Ax 

Interaction eax 

16 

particle-wall 0.075 

16 

particle-particle 0.025 

32 

particle-wall 0.05 

32 

particle-particle 0.025 
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FIG. 5. (Color online) Lubrication corrections for the cases of normal particle-wall (a) and particle-particle 
(b) interactions, respectively compared against the analytical solution of Brenner 


Results shown for two different resolutions, Dp/lS.x = 16 and 32. 


^1 and Cooley and ONeill 


below a threshold Ea so that A(e < e^) = A(eo-). This threshold value is related to the typical size 
of the asperities and was fixed to £„ = 0.001 for particle-wall interactions. We keep the Stokes 
amplification factor saturated until the surfaces overlap; then the collision force takes over. Hence, 
the force acting on the particle is corrected by given by: 


^Fiub 
QtT /iRpUij^Yi 


A(e) - A(eAx), <e < eax 
< A(e^) - A(eAx), 0 < e < 


(34) 


0, otherwise. 

For particle-wall collisions, the normal fluid-induced forces are set to zero for overlaps larger 
than the overlap due to the particle’s submerged weight, = \pp — pf\gVp/kn, in order to avoid 
artificial dissipation due to the stretching of the collision time of the contact model. This procedure 
is not extended to particle-particle interactions, as it can canse significant artificial increase in the 
particles’ acceleration for colliding particle pairs dne to a sndden decrease in drag force. 


B. Piezoviscous effects 


Joseph and Hunt 


23l | performed experiments on wet, oblique collisions of spheres onto planar 


surfaces. They showed that the coefficient of sliding friction decreased by one order of magnitnde 
when their smooth steel spheres collide, whereas it remained of the same order of magnitude 
(~ 15% higher) for the case of rough glass spheres. They suggested that this abrupt decrease of 
the friction coefficient for smooth spheres was due to the fact that a characteristic piezoviscous 
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lenghtscale 


36l |. hpy, was larger than the average size of the asperities and therefore ’contact’ occurs 


through the fluid, which is behaving like an elastic-solid. Also, the slight increase in the coefficient 
of sliding friction for rough spheres was explained by the fact that the fluid introduces an extra 
resistance when the asperities have relative motion in the tangential direction. They developed 
a model capable of predicting the coefficient of sliding friction of smooth spheres colliding onto 
planar surfaces in a viscous liquid from elastohydrodynamic lubrication theory. 

Hence, for the case in which piezoviscous effects are important, it does not suffice to use input 
parameters from dry collisions and lubrication corrections for obtaining a physically realistic result: 
the coefficient of sliding friction measured in a wet collision experiment, or predicted by the model 
developed in fic,wet, should be used. 


IV. NUMERICAL IMPLEMENTATION 

The fluid phase is governed by the Navier-Stokes equations for an incompressible Newtonian 
fluid: 

V • u = 0, (35a) 

O -I 

-I- V • uu =-Vp -|- i/f V^u, (35b) 

dt pf 

where u is the fluid velocity, Uf = p/pf the kinematic viscosity of the fluid and p the pressure. 

The translational and rotational motion of solid particles is described by the Newton-Euler 
equations for rigid body motion. For a spherical particle they read, 

du. f 

Pp^p~^ = f T ■ ndA +{pp - pf)Vpg + Fc, (36a) 

dt Jdv 

Ip— 7 ^ = (f r X (r • n) dA + Tc- (36b) 

dt Jdv 

The left-hand-side of (I36ap is the temporal variation of linear momentum of the particle: Uc is the 
centroid velocity, pp the particle mass density, and Vp the particle volume given by (4/3)7ri?p for a 
spherical particle with radius Rp. 

The first term of the right-hand-side of (|36ap is the net force resulting from the distribution of 
fluid stress, r = —pl+p (Vu -|- Vu'^) at the particle surface, dV, projected to the outward-pointing 
unit normal to dV, n. The second term is the buoyancy force due to a difference between the fluid 
and particle densities in the presence of a gravitational field with acceleration g. Fc represents 
other external forces acting on the particle such as, e.g., collision forces. 

The left-hand-side of (I36bp is the temporal variation of angular momentum of the particle, where 
Uc is its angular velocity and Ip the moment of inertia of a solid sphere, given by {2/5)ppVpRp. 
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Due to spherical symmetry, only two non-trivial terms balance the left-hand-side: the flow-induced 
torques and the external torques (e.g., a collision torque) respectively the first and second terms 
in the right-hand-side of (|36b(l . r = x — Xc is the position vector relative to the particle centroid 
X = Xc- Tc is an external torque that acts on the particle whenever there is a contact force with 
a tangential component. 

The equations (I35ap , (I35bp and (I36al) , (I36bp form a set of equations coupled through no-slip and 
no-penetration (ns/np) boundary conditions at the particle/fluid interface. Hence, the velocity at 
the particle surface, 

Up = Uc cJc X r, (37) 

is required to match the local fluid velocity: 

u = Up (x) V X € dV. (38) 


The governing equations for the fluid phase are integrated in time with an explicit low-storage 
three-step Runge-Kutta method for all terms except the pressure gradient in the Navier-Stokes 
equations; for the latter the Crank-Nicolson scheme is used. The equations are discretized in space 
on a uniform, staggered Cartesian grid with the finite-volume method in which spatial derivatives 
are estimated with the central-differencing scheme. To enforce ns/np conditions at the particle’s 
surface we use the second-order accurate IBM developed in j^. The advantage of an IBM is 
that the governing equations are solved on a spatially continuous grid without any holes, which 


enables the use of an efficient, FFT-based, direct solver for the pressure Poisson ec 
restrictions for the computational time step have been derived by Wesseling 


nation. Stability 


37l |. For a uniform 


Eulerian grid with Ax = A?/ = Az and the central-differencing scheme, a sufficient criterion for 
von Neumann stability is given by: 

At 


Cou = 


1.65 y/SAx 

12 ’ eTTh 


min 


< 1 . 


(39) 


The governing equations for the solid particles are advanced in time with the same Runge-Kutta 
scheme as used for the fluid phase, except for collision forces/torques and tangential displacement; 
these terms are integrated with a second-order Crank-Nicolson (CN2) scheme that has proven to 
return a stable and accurate integration. This scheme requires the contact force at the next time 
level, q, which depends on the values of the particle position and velocity at the same level (Eqs. (|15l) 
and (1221) ). We therefore compute the contact force iteratively as a function of the particle position 
and velocity at q until the new particle position converges. The particles’ position and velocity 
are initialized {k = 0) with the values of the previous time level q — 1- The advancement follows 
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directly the integration of the Navier-Stokes equations, within the RK3 time advancement loop 
with a time step Atp which is allowed to be smaller than the time step of the Navier-Stokes solver 
At to ensure that the contact forces and lubrication force corrections are accurately integrated. 
The forces induced by the IBM are fixed in time while the sub-integrations are performed. For 
sub-stepping ratios r^t = At/Atp ranging from 1 to 0(100), the extra overhead introduced by the 
sub-stepping is negligible. The scheme is illustrated below, 
fc = 0 
do 

for all particles j in contact with particle i do 

compute 61^% and = R • S'}-/ + + S ' uf”"/) 

compute FFj; and 
update and 

end for 


= u: 


9-1 


(particle-fluid coupling terms 
Atl 


m 


At} Fg’fe -fFg-i 

2 PpFp 


^ 1 / = xr 1 (u«’'=+ur 1) 


= a; 


9-1 


-I- (particle-fluid coupling terms 


Atl -f 


Ppip 


err,t„^ = x^’ - x. 


9,fe _ ^9.^-11 


(40) 

(41) 

(42) 

(43) 


k = k + 1 

while < erruer,max 

Atp = (ur -f (3r)Atp varies accordin g to the duration of the Runge-Kutta sub steps and the 

ai = 32/60, /?i = 0, ^2 = 25/60, /?2 = -17/60, 


coefficients can be found in Wesseling 


03 = 45/60, = —25/60. The lubrication corrections are integrated with the same scheme 

as the collision force. From this CN2 scheme, we expect second-order accuracy for the linear 
momentum of the particle and consequently, third-order accuracy for the integration of the particle 
velocity. We verified the accuracy of the method by reproducing, in simulations of dry collisions, 
the coefficients of restitution and et,d) that are used as an input in the collision model (not 
shown). For the simulations of the present work 1 iteration sufficed for obtaining a small iterative 
error: errj^^j. > 10“®Aa:. 

Unless otherwise stated, the particles are resolved with Dp/Ax = 16 and a sub-stepping ratio 
of r^t = 50, the collision time set to Tn = 8At and the time step set by Cou = 0.5. 
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V. RESULTS FROM COLLISION SIMULATIONS 


A. Bouncing motion of a solid sphere colliding onto a planar surface in a viscous liquid 


We simulated the bouncing motion of a solid sphere immersed in a viscous liquid and colliding 
under gravity onto a planar surface. The trajectory of the point of the particle closest to the 
surface and time evolution of its velocity are compared to the experimental data of Gondret et al. 
Q]. This experiment is a useful benchmark for confirming that the lubrication corrections and 
collision model return a realistic bouncing velocity, and that the collision is represented in good 
approximation as an instantaneous event in the particle motion. Furthermore, there is no need 
for specifying impact and rebound velocities, which definitions vary significantly in literature Q- 
Note that small differences in rebound velocity are amplified after its temporal integration, and 
therefore more noticeable in the particle trajectory. Hence, a good agreement with this experiment 
gives a fine indication that the approach used to resolve a head-on wet collision is adequate. 

The simulations were carried out in a domain corresponding to a closed container with dimen¬ 
sions Lx/Dp X Ly/Dp X Lz/Dp = 12 x 30 x 12. The particle is initially placed at y/Dp = Ly — l.^Rp, 
centered in Lxj^ and Lzl2. The motion is driven by a downward-pointing gravitational accelera¬ 
tion of <7 = 9.81m/s^. The time step was fixed to the maximum allowed by the stability criterion at 
the maximum particle velocity (i.e., at impact), multiplied by Cou to ensure a stable and accurate 
temporal integration. The physical and computational parameters are listed in Table ITIl 


Figure 6(a) presents the results for the trajectory and time evolution of velocity of a steel sphere 
colliding onto a glass wall immersed in silicon oil RVIO, corresponding to Case St,i = 152 of Table 
m The model is able to accurately reproduce this case. Moreover, the large discrepancy for the 




(a) (b) 

FIG. 6. (Color online) Trajectory (a) and time evolution of the particle velocity (b) in the bouncing motion 
of a steel sphere colliding onto a planar surface in silicon oil RVIO. 
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numerical solution in the absence of lubrication model illustrates the importance of including it. 
Note that at each impact the particle has a lower Stokes number: St^pstb = 152, St„_ 2 ndb = 81, 
St,i^3rdb — 23 and St,^^4thb — 10. 

Figure 7(a) compares our simulations to the experimental data of Gondret et al. [91] of the 


first bounce of steel spheres colliding onto planar surfaces in silicon oil at different impact Stokes 
numbers. In the cases for which Ly was not sufficiently large for the particle to reach its terminal 
velocity before colliding with the wall, we imposed an initial velocity to the particle to ease the 
convergence of the velocity to its terminal value. For extreme cases of a highly inertial St„ = 742 
and highly viscous St^ = 29 flow the resolution was increased to Dp/Ax = 32. 


TABLE II. Properties of the fluids and solid spheres used in the experiment of Gondret et al. |9| and 
computational parameters of the numerical simulations. 


Case 

Dp [mm] 

Pp [kg/m^j 

^n,d 

P [cP] pf [kg/m^] 

Dp/ Ax 

Cou 

rAt 

N 

St„ = 742 

5 

7800 

0.97 

5 

920 

32 

0.5 

50 

8 

St„ = 152 

3 

7800 

0.97 

10 

935 

16 

0.2 

50 

8 

St„ = 100 

4 

7800 

0.97 

20 

953 

16 

0.2 

50 

8 

St„ = 29 

6 

7800 

0.97 

100 

965 

32 

0.5 

50 

8 




t -tc{s) t - tc (s) 


(a) 


(b) 


FIG. 7. (Color online) Trajectories obtained from simulations of particles colliding onto a planar surface in 
silicon oil, for different impact Stokes numbers with (a) and without (b) closure for lubrication interactions. 
Experimental data from Gondret et al. Q. 


As expected and shown in Figure [7(b) [ the deviation from the experimental data for the simula¬ 
tions without lubrication closure is more significant for smaller Stokes numbers due to the increasing 
importance of viscous effects. The simulations show a good agreement with the experimental data 
for this wide range of Stokes numbers. 
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Sensitivity of the results to the collision time and, time step and sub-stepping 


In the following we explore the sensitivity of the model to the computational parameters that 
govern the collision time and temporal integration of the fluid and particle motion. These param¬ 
eters are prescribed collision time, amount of sub-stepping and time step of the overall numerical 


algorithm. Let us consider the trajectory of Figure 6(a) as the reference case for this sensitivity 
analysis, with focus on the first bounce (the subsequent will be influenced by how realistically the 
first is reproduced). We performed a set of simulations with parameters shown in Table Hill 


TABLE III. Computational parameters used for the sensitivity study. 


Case 

Cou 

VAt 

N Sffff/Ax 

(%) Notes 

REF 

0.2 

50 

8 

33.6 

Reference case 

SAl 

0.6 

50 

5 

63.0 

Larger At, smaller N 

SA2 

0.2 

1 

8 

33.6 

No sub-stepping 

SA3 

0.2 

50 

1 

4.16 

Small N 

SA4 

0.025 

50 

8 

4.16 

Smallest At, same T„ as SA3 

SA5 

0.1 

50 

16 

33.6 

At between REF and SA5, same T„ as REF 

SA6 

0.025 

50 

64 

33.6 

Same At as SA4, same T„ as REF 


Figure 8(a) presents the outcome of this set of simulations. The trajectory corresponding to case 
SAl compares well with the one of REF, which shows that a collision which takes 5 Navier-Stokes 
can still be realistically reproduced. 

The trajectories of cases SA2 and REF cannot be distinguished; this shows that sub-stepping 
is not required to better resolve the collision and lubrication force corrections in this case because 
the time step of the Navier-Stokes solver is sufficiently small. 

In case SA3 the collision time takes exactly one time step of the Navier-Stokes solver, which has 
the same value that the one of REF. Although the collision force and lubrication corrections are 
resolved due to the sub-stepping, the trajectory obtained from this simulation differs significantly 
from the reference case. This is mostly a consequence of an over-estimation of the drag force from 
the IBM when the surrounding fluid does not adapt itself gradually to the abrupt change in particle 
velocity due to a collision, as illustrated hereafter. 


Decreasing the time step of SA3 while keeping the stiffness fixed (SA4) allows the fluid to adapt 
itself to the changes in particle velocity. However, the simulation also over-estimates the drag force 
acting on the particle. We further show with cases SA5 and SAG that the over-estimation of the 
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drag force is not consequence of an inconsistency problem, because the simulations, for the same 
particle stiffness, converge monotonically to SAG with decreasing time step. 

The discrepancy of the solution for the stiff particle of case SA4 is caused by a loss of conservation 
properties of the interpolation kernel used by the IBM when its stencil, for a certain Lagrangian 


forcing point, overlaps with the one of another particle or with a solid wall 


381. This issue 


becomes signihcant for considerably high particle stiffnesses, where more problematic forcing points 
continue to perform interpolation/spreading operations in a inconsistent manner throughout the 


entire collision time. Figure 8(b) shows simulations for cases SA3* and SA4* with the same 
parameters as the ones of SA3 and SA4, but excluding from the forcing scheme Lagrangian forcing 


points with a distance to the wall smaller than Ax (procedure similar to what is suggested in 381] ) ■ 


Indeed, simulation SA3* still yields an over-estimated drag force, whereas SA4* yields the realistic 
bouncing trajectory with a difference in the peak of the trajectory of 2.5% from REF. 

This illustrates that the realistic bouncing trajectory can only be obtained if the surrounding 
fluid is allowed to adapt itself to the changes in particle velocity. Hence, we decided to ensure that 
the fluid phase adapts itself to the changes in particle velocity by avoiding excessively high values 
of particle stiffness. Note that for the reference case the maximum overlap is already significantly 
small, about one third of a grid cell. 




FIG. 8. (Color online) Sensitivity analysis to the time step, sub-stepping and stretching of the collision time 
(a) and outcome of cases SA3 and SA4 when problematic Lagrangian forcing points are excluded from the 
IBM forcing scheme (b). Computational parameters in Table lllll 


B. Wet head-on collisions 


The previous validation gives a fine indication of the realism of the approach used to simulate a 
wet head-on collision. On contrary, the experimental curves of enl&n,d = /(St) - benchmark often 
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used to validate these models - depend on the definition of impact and rebound velocities that are 
used to compute Cn- If, for instance, we define Uin,n as the terminal velocity, and Uout,n as the 
maximum velocity after impact, for the case St^ = 152 of Table HIl we obtain = 0.85; considerably 
different from the experimentally measured value of 0.78. To circumvent this problem one can define 
impact and rebound velocities which agree with the frame rates used in the measurements [l^ . We 
therefore use the impact velocity and rebound velocities at the instants t — tc = T/~^, respectively; 
where / is a frequency related to the temporal resolution of the experiment. 


Particle-wall collisions 


We simulated particle-wall collisions in a viscous liquid for several values of St^ and cornpared 
the resulting normal coefficients of restitution e„ to the experimental data of Joseph et al. 33l |. 

The computational domain has dimensions of L^/Dp x Ly/Dp x L^/Dp = 12 x 24 x 12. Similarly 
to the previous cases, the particles are placed at a distance y/Dp = Ly — 1.5Rp and their motion 
driven by gravity. We simulated steel spheres colliding onto a planar surface in silicon oil RV20 
(physical parameters are listed in Table HH) . The Stokes number was varied by varying the particle 
diameters from 1.5 mm to 10 mm. We used a value of / = 500 Hz, which complies with the 
frequency of image acquisition of the experiment. Figure [9] shows the results. 
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FIG. 9. (Color online) Normal, wet coefficients of restitution for particle-wall collisions. The experimental 
data were normalized with the value e„ d = 0.97 measured in the reference. 


The numerical simulations agree with the experiments for the entire range of impact Stokes 
numbers. This agreement is expected after the finer validation of the previous section, and careful 
definition of impact and rebound velocities. 
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Particle-particle collisions 


For inter-particle collisions we reproduced the pendulum experiment of Yang and Hunt 


2 ^ by 


colliding a moving projectile particle with a steady target particle. Spheres of the same size and 
material were centered in a computational box with dimensions Lx/Dp x Ly/Dp x Lz/Dp = 6x12x6 


and separated in the y-direction by a distance of 4 Dp. Similarly to Simeonov and Calantoni 


391. we 


force an acceleration g to the projectile particle to mimic the release mechanism of the experiment. 

The physical parameters are comparable to the experiments of head-on collisions of steel spheres 
in aqueous solutions of glycerol: Cn h = 0.97, pp = 7780kg/m^, Dp = 12.7mm, p = 45 cP, 
Pf = 1125kg/m^. Yang and Hunt 25| defined the rebound and impact velocities at instants 
corresponding to a value of / = 100 Hz. 

The binary impact Stokes number, defined as Stij^n = {^/^)PpUij^nDp/p for two equal spheres of 
the same material, was changed by varying the projectile particle’s acceleration from g = 0.02 x 9.81 
to 10 X 9.81 m/s^. We used a value of Eo- = 10“® to resolve the lubrication interaction in the thin 
gap-width between these smooth particles. This value agrees with the order magnitude of the 
size of the asperities (0(0.1) — 0(0.01) pm 25|). These small values together with the fact that 


the target particle is freely mobile (numerical solution more sensitive to errors when compared 
to a collision with a wall or a fixed particle) make this benchmark a valuable test for the overall 
methodology. Resolving the lubrication layer of the interacting particles at such a small scale 
required a time step dictated by Cou = 0.1 for a resolution of Dp/Ax = 16, and a sub-stepping 
ratio of rAt = 1000. For values of St,i higher than 0(100), the resolution required to describe the 
dynamics of the intervening film is higher. Hence Dp/Ax was increased to 32, with a time step 
dictated by Cou = 0.5. 


Figure 10(a) presents the trajectories of the particles’ contact points (results of the numerical 
simulations shifted vertically for clarity). For very small impact Stokes numbers, the momentum 
transferred to the target particle is not sufficiently high for it to overtake the viscous drag and 


travel independently. Yang and Hunt 


251] observed that this is the case for Stij^„ < 10, where the 


particles tend to move as a pair with constant separation distance. This is shown in Figure 10(a 


for cases St = 11.8 (measured experimentally) and St = 12.7 obtained from a numerical simulation. 
The good agreement between the numerical simulation and the experiment is a strong indicator 
of the success of the overall method to resemble this viscous limit. In particular, it gives a finer 
assessment of the realism of the lubrication closure. Furthermore, the simulations with values of 
binary impact Stokes number considerably larger than 10, St = 21.5 and St = 34.3 do not show 
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this trend, which is consistent with the experimental observations. 

Finally, Figure [10] compares the computed effective binary coefficient of restitution from the 
numerical simulations to the experiments. The necessity of increasing the spatial resolution of 
the simulation for a binary impact Stokes number of St = 135 is also illustrated by showing the 
outcome of this case with both resolutions. Increasing the resolution becomes more important in 
this case than in particle-wall interactions due to the requirement of an accurate description of the 
interacting dynamics of the two particles through short-range hydrodynamic interactions. 




FIG. 10. (Color online) Trajectories of the particles’ contact points(results of the numerical simulations were 
shifted vertically for clarity). The solid line was extracted from [23 (a). Wet coefficients of restitution for 
particle-particle collisions (b). The experimental data were normalized with the value en,d = 0.97 measured 
in the reference. 


The agreement with the experimental data further supports the validity of our approach. We 
should note that extra computational overhead (Con = 0.1 for Dp/Ax = 16) was required for 
reproducing these results, when compared to particle-wall collisions. 


C. Oblique collisions 


Finally, we validated our model for oblique particle-wall co. 
liquids. We use the experimental data of Joseph and Hunt 


lisions in a dry system and in viscous 


23[ of oblique particle-wall collisions 


in air and aqueous solutions of glycerol. The collisional properties parameters of the particles 
agree with their experiments and are described together with the other physical parameters of the 
simulations in Table nYi The computational domain and particle’s initial position is the same of 
the previous simulations of particle-wall collisions. The particle motion is driven by an imposed 
acceleration with direction eg = —sm{4>in)ey — cos{4>in)ez, to yield the desired incidence angle. 
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The magnitude of the particle acceleration was set to g = 10 x 9.81 m/s^ to ensure that the 
glass spheres collide with an impact Stokes number of 0(1000), comparable to the experimentally 
measured values. The results for immersed collisions of steel spheres show little sensitivity to the 
choice of the value of the acceleration due to the small value of the coefficient of sliding friction. 

TABLE IV. Physical and computational parameters for the simulations of oblique particle-wall collisions. 

JVIaterial Dp Gn,d fdc,wet Pp Pf P 

steel 2.5 mm 0.97 0.34 0.11 0.02 7800 kg/m^ 998 kg/m^ 1 cP 

glass 2.5 mm 0.97 0.39 0.10 0.15 2540 kg/m^ 998 kg/m^ 1 cP 


Figure [TT] shows a comparison between the normalized incidence and rebound angles obtained 
from oblique collisions between steel and glass spheres. 
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FIG. 11. (Color online) Results of for oblique collision simulations in a dry system (a), and in a viscous 
liquids (b). Experimental data of Joseph and Hunt 231. 


The simulations agree well with the experimental data for the entire range of incidence angles. 
This is an expected consequence of the fact that the model uses the macroscopic properties of these 
collisions as input parameters. 


VI. CONCLUSIONS AND OUTLOOK 

We presented and validated a collision model for fully-resolved 4-way coupled simulations of 
flows laden with finite-size solid particles. There are three types of particle-particle or particle-wall 
interactions that must be reproduced in such simulations: (1) long-range hydrodynamic interac¬ 
tions; (2) short-range hydrodynamic interactions; and (3) solid-solid contact. 

The long-range hydrodynamic interactions are computed by a Navier-Stokes solver where we 
used an IBM for an efficient representation for the particles. Other approaches that require a closure 
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for small inter-particle/particle-wall distances (e.g., Lagrangian-multiplier or Lattice-Boltzmann 
methods) could have also been used. 

Short-range hydrodynamic interactions are also partly resolved by the IBM. However, the dis¬ 
crete nature of these numerical methods together with the necessity of a computationally efficient 
implementation typically require a closure model for lubrication interactions. For the cases ad¬ 
dressed here, the only lubrication interaction that requires modeling is the squeezing of fluid through 
the thin gap between two approaching particles or a particle approaching a wall. To achieve this 
we used a two parameter model: for normalized gap-widths smaller than a value gax we introduce 
a correction based on asymptotic expansions of analytical solutions of particle-particle/-wall inter¬ 
actions in the Stokes regime. This value is obtained by determining the gap-width for which our 
numerical method is unable to reproduce the lubrication interaction. The second parameter, 
accounts for roughness effects for even smaller gap-widths. 


Finally, solid-solid contact is modeled through a linear soft-sphere collision model capable of 
stretching the collision time, to avoid computational overhead in the calculation of the collision 
force. The model constants are analytically related to the three input parameters of the model 
described by Walton 28|, which are widely reported in the literature. The model can be extended 


to accommodate more complex mechanics such as adhesion or plasticity for the normal force, or 
static and dynamic friction for the tangential force. However, these features are in general not 
required in 4-way coupled simulations of flows with finite-size particles at small/moderate solid 
volume fractions. 


We validated our methodology against several benchmark experiments and the results show a 
good quantitative agreement. The simulations of the bouncing trajectory of a spherical particle 


colliding onto a planar surface ^ show that the lubrication force corrections, combined with the 
collision model are sufficient for reproducing a realistic bouncing velocity. Subsequently, we suc¬ 
cessfully reproduced experimental data for the normal coefficient of restitution as a function of the 


impact Stokes numbers for head-on particle-wall 


33l | and particle-particle collisions 2^. Finally, 


our simulations of oblique particle-wall collisions in dry and wet systems agree quantitatively with 
the experimental data of Joseph and Hunt 2^ for the entire range of incidence angles. We reserve 


further validations of the overall model for flows with many particles where both lubrication and 
friction play an important role for a future publication. 


The physical realism and computational efficiency of our method allows for massive fully- 


resolved simulations of particle-laden flows with 4-way coupling. 
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